#Selecting the genetic instruments for traits of interest

library(tidyverse)#v2.0.0
library(bigreadr)#v0.2.5
library(LDlinkR)#v1.3.0

#example with GSCAN Smoking Initation: https://conservancy.umn.edu/bitstream/handle/11299/201564/SmokingInitiation.txt.gz?sequence=34&isAllowed=y

#format GWAS summary statistics with following fields: CHR,BP,SNP,A1,A2,N,SE,P,OR|BETA,MAF,INFO
sumstats<-fread2('./Data/SmokingInitiation.txt.gz')# this file can be downloaded from the link above

sumstats<- sumstats %>% dplyr::select(CHROM,POS,RSID,ALT,REF,EFFECTIVE_N,SE,PVALUE,BETA,AF) %>% mutate(INFO=1)

colnames(sumstats)=c('CHR','BP','SNP','A1','A2','N','SE','P','BETA','MAF','INFO')

#keep biallelic SNPs with MAF>1%
Multi_allelic<- sumstats %>% group_by(SNP) %>% summarize(count_a0=length(unique(A2)),count_a1=length(unique(A1))) %>%
  filter(count_a1>1 | count_a0>1) %>% select(SNP) %>% mutate(rm_reason='Multi-allelic')
Indels<- sumstats %>% filter(nchar(A1)>1 | nchar(A2)>1) %>% select(SNP) %>% mutate(rm_reason='INDEL')
Rare<- sumstats %>% filter(MAF<0.01 | MAF>0.99) %>% select(SNP) %>% mutate(rm_reason='RareSNP')
sumstats<-filter(sumstats, ! SNP %in% c(Multi_allelic$SNP,Indels$SNP,Rare$SNP,'.'))#filtered GWAS sumstasts

#1.Run clumping on at r2=0.01 with 10MB LD windown in European 1000 Genomes plink files
#example with 1000 Genomes genotype File: ./Data/QC_1000G_P3* (VCF files downloaded from https://www.internationalgenome.org/category/genotypes/)
#plink (v1.9b) command to convert vcf to bed files (plink --vcf vcf_name --make-bed --out QC_1000G_P3)
#bash script have 4 args: output_directory path, trait, if binary, stats (if BETA or OR)  

system("sbash ./PRS_1000G_byTraits.sh ./Data/ trait true BETA")

#2.Import clumped rsid list in 1000G and check if in ccRCC genotype file 

list_vars<-read_table2("./Data/SmokingInitiation_1000G_0.01.snp")#list clumped variants
list_vars<-filter(list_vars,P<=5e-08)# keep Genome-wide hits
list_ccRCC_vars<-fread2('./Data/ccRCC_Genotype.txt')#list variants found in ccRCC genotype plink file
list_ccRCC_vars %>% filter(V2 %in% sumstats$SNP)->list_ccRCC_vars
notfound<-list_vars %>% filter(!SNP %in% list_ccRCC_vars$V2)

#3.LD proxy for missing variants in ccRCC genotype file 
#need to register to get token to use this function

ldproxy=data.frame()
for(i in notfound$SNP){
  skip_to_next <- FALSE
  temp_data<-LDproxy(i, pop = "EUR", r2d = "r2", token = '', file = FALSE)
  tryCatch(temp_data %>% select(RS_Number,R2) %>% 
             filter(RS_Number %in% list_ccRCC_vars$V2 & R2>=0.8) %>% 
             arrange(R2) %>% slice(1) %>% 
             mutate(query=i, ProxyFound='Yes')-> temp_data, error = function(e) { skip_to_next <<- TRUE})
  if(skip_to_next) { next }
  colnames(temp_data)=c('proxy_SNP','R2','SNP','ProxyFound')
  ldproxy=rbind(ldproxy,temp_data)
}

proxyfound=nrow(ldproxy)#n proxy found
proxyNotFound=data.frame(proxy_SNP=NA,R2=NA,SNP=setdiff(notfound$SNP,ldproxy$SNP),ProxyFound='No')# proxies not found
ldproxy=rbind(ldproxy,proxyNotFound)

#4.Update basefile with proxies
sumstats %>% filter(SNP %in% list_vars$SNP | SNP %in% ldproxy$proxy_SNP) %>%
  mutate(Is_Proxy=ifelse(SNP %in% ldproxy$proxy_SNP, 'Yes','No'))->basefile

#5.Variance of trait explained by genetic 
basefile<- mutate(basefile, VE=((BETA^2*2*MAF*(1-MAF))/(BETA^2*2*MAF*(1-MAF)+SE^2*2*N*MAF*(1-MAF))))
sumstats %>% filter(SNP %in% proxyNotFound$SNP)->basefile_missing
basefile_missing<- mutate(basefile_missing, VE=((BETA^2*2*MAF*(1-MAF))/(BETA^2*2*MAF*(1-MAF)+SE^2*2*N*MAF*(1-MAF))))

#6.Create table summarising variants per trait
table_variants=data.frame(GWAS_trait='trait',
                          N_var_1000G=nrow(list_vars),
                          Found_Needed=paste0(proxyfound,"/",nrow(notfound)),
                          Final_var=nrow(basefile),
                          VE=sum(basefile$VE)*100,
                          VE_miss=sum(basefile_missing$VE)*100)

#7.Save files
write_csv(table_variants,paste0('./Data/SmokingInitiation_summary_SNPs.csv'))
write_tsv(basefile,paste0('./Data/SmokingInitiation_summary_SNPs.txt'))#Example_Basefile_SmokingInitiation.txt
write_tsv(basefile_missing,paste0('./Data/SmokingInitiation_summary_SNPs_missing.txt'))

#8. Calculate Polygenic Risk Score (PRS) in ccRCC dataset
#bash script have 4 args: output_directory path, trait, if binary, stats (if BETA or OR)  
system("sbash ./PRS_ccRCC_byTraits.sh ./ trait true BETA")

